
# Preliminaries -----------------------------------------------------------

rm(list=ls())

# Set directories ---------------------------------------------------------

# get user's m1
home_directory = "[your home directory]" # change this to your home directory
master_scripts_directory <- file.path(home_directory, "/Figures/Non_IRF_Figures") # where the program sources DFA functions for figures
figures_directory <- file.path(home_directory, "/Figures/Non_IRF_Figures/output")

# Load in data ------------------------------------------------------------
library(RColorBrewer)
library(latex2exp)
library(readxl)
library(tidyverse)
source(file.path(master_scripts_directory, "functions.R"))
source(file.path(master_scripts_directory, "dfa_plot_funcs.R"))

### This script gathers the data needed for all of the charts 
source(file.path(master_scripts_directory, "get_businesscycle_data.R"))

buscycle_legtext <- c("Top .1%", "99 - 99.9%", "90 - 99%", "70 - 90%", "50 - 70%", "Bottom 50%")
buscycle_colors <- rev(dfa_colors[1:6])

setcmpar()
style$legend.cex <- 1
style$freq.label.cex <- 1

# Appendix figures --------------------------------------------------------


# # Figure 15: Alternative error assumptions --------------------------------
# 
# lit_shares <- read_dfa(.item = "shares", .group = "networth", .agg_level = "bus_cycle2", .method = "litterman") %>%
#   mutate_if(is.numeric, ~.*100) %>%
#   select(date, cat, networth) %>%
#   rename("litterman" = networth)
# 
# cl_shares <- read_dfa(.item = "shares", .group = "networth", .agg_level = "bus_cycle2", .method = "chowlin1") %>%
#   mutate_if(is.numeric, ~.*100) %>%
#   select(date, cat, networth)  %>%
#   rename("chowlin" = networth)
# 
# shares_comp <- Reduce(left_join, list(select(dfa_shares, date, cat, networth) %>% rename("fernandez" = networth), lit_shares, cl_shares)) %>%
#   pivot_wider(id_cols = date, names_from = cat, values_from = c(fernandez, litterman, chowlin)) %>%
#   mutate(date = as.Date(date, frac=1))
# 
# shares_comp_top <- shares_comp %>%
#   select(date, ends_with("TopPt1"), ends_with("pct99to999"), ends_with("pct90to99"))
# 
# postscript(file = file.path(figures_directory, "shares_comp_top.eps"), width = 4, height = 4)
# par(mar = c(2,1,2,2))
# pinfo <- rplot.line(
#   plotdata = shares_comp_top,
#   Col = rep(c("black", "red", "blue"), 3),
#   Y2lab = "Percent of aggregate",
#   Enddatelab = FALSE,
#   Slabvec = c("Top .1%", "99 - 99.99%", "90 - 99%"),
#   Slabloc.x = 2021,
#   Slabloc.y = c(8, 20, 40),
#   Slabadj = c(1, 0),
#   Slabcol = "black"
# )
# dev.off()
# 
# shares_comp_bottom <- shares_comp %>%
#   select(date, ends_with("pct70to90"), ends_with("pct50to70"), ends_with("Bottom50"))
# 
# postscript(file = file.path(figures_directory, "shares_comp_bottom.eps"), width = 4, height = 4)
# par(mar = c(2,1,2,2))
# pinfo <- rplot.line(
#   plotdata = shares_comp_bottom,
#   Col = rep(c("black", "red", "blue"), 3),
#   Y2lab = "Percent of aggregate",
#   Enddatelab = FALSE,
#   Slabvec = c("70 - 90%", "50 - 70%", "Bottom 50%"),
#   Slabloc.x = 2021,
#   Slabloc.y = c(23, 9, 3),
#   Slabadj = c(1, 0),
#   Slabcol = "black"
# )
# dev.off()
# 
# postscript(file = file.path(figures_directory, "shares_comp_legend.eps"), width = 6, height = .6)
# par(mar=c(0,0,0,0))
# plot(NULL ,xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
# legend("top", legend = c("Fernandez", "Litterman", "Chow-Lin"), cex=1.25,
#        bty='n', col = c("black", "red", "blue"), horiz = TRUE,
#        lty = 1, lwd = 2)
# dev.off()


# Figure 16a: standard errors of estimates ------------------------------------------------

postscript(file = file.path(figures_directory, "shares_se_top.eps"), width = 4, height = 4)
par(mar = c(2,1,2,2))
pinfo <- rplot.line(
  plotdata = fernandez_nw_shares_withse_top,
  Col = c(rep(buscycle_colors[6], 3),  rep(buscycle_colors[5], 3),  rep(buscycle_colors[4], 3)),
  Lty = c(rep(c("solid", "dashed", "dashed"), 3)),
  Y2lab = "Percent of aggregate",
  NBER = TRUE,
  shading.border = FALSE,
  Enddatelab = FALSE,
  Y2lim = c(-2,45),
  Y2at = seq(0, 45, 5),
  Zeroline = TRUE
)

legend(x = "top", 
       legend = c(buscycle_legtext[1:3]),
       col = c(buscycle_colors[6:4]),
       x.intersp = .3,
       bty = "n",
       cex = .65,
       lwd = 1.5,
       ncol = 3)
dev.off()

postscript(file = file.path(figures_directory, "shares_se_bottom.eps"), width = 4, height = 4)
par(mar = c(2,1,2,2))
pinfo <- rplot.line(
  plotdata = fernandez_nw_shares_withse_bottom,
  Col = c(rep(buscycle_colors[3], 3),  rep(buscycle_colors[2], 3),  rep(buscycle_colors[1], 3)),
  Lty = c(rep(c("solid", "dashed", "dashed"), 3)),
  Y2lab = "Percent of aggregate",
  NBER = TRUE,
  shading.border = FALSE,
  Enddatelab = FALSE,
  Y2lim = c(-2,45),
  Y2at = seq(0, 45, 5),
  Zeroline = TRUE
)

legend(x = "top", 
       legend = c(buscycle_legtext[4:6]),
       col = c(buscycle_colors[3:1]),
       x.intersp = .3,
       bty = "n",
       cex = .65,
       lwd = 1.5,
       ncol = 3)
dev.off()


# Figure 16b: Kalman results -------------------------------------------------------------

# Kalman with standard errors, layered over baseline shares 

postscript(file = file.path(figures_directory, "cred_comp_top.eps"), width = 4, height = 4)
par(mar = c(2,1,2,2))
pinfo <- rplot.line(
  plotdata = kal_nw_shares_comp_withse_top,
  Col = c(rep(buscycle_colors[6], 3), "black", rep(buscycle_colors[5], 3), "black", rep(buscycle_colors[4], 3), "black"),
  Lty = c(rep(c("solid", "dashed", "dashed", "solid"), 3)),
  Y2lab = "Percent of aggregate",
  NBER = TRUE,
  shading.border = FALSE,
  Enddatelab = FALSE,
  Y2lim = c(-2,45),
  Y2at = seq(0, 45, 5),
  Zeroline = TRUE
)

legend(x = "top", 
       legend = c(buscycle_legtext[1:3], "Baseline"),
       col = c(buscycle_colors[6:4], "black"),
       x.intersp = .3,
       bty = "n",
       cex = .65,
       lwd = 1.5,
       ncol = 2)
dev.off()

postscript(file = file.path(figures_directory, "cred_comp_bottom.eps"), width = 4, height = 4)
par(mar = c(2,1,2,2))
pinfo <- rplot.line(
  plotdata = kal_nw_shares_comp_withse_bottom,
  Col = c(rep(buscycle_colors[3], 3), "black", rep(buscycle_colors[2], 3), "black", rep(buscycle_colors[1], 3), "black"),
  Lty = c(rep(c("solid", "dashed", "dashed", "solid"), 3)),
  Y2lab = "Percent of aggregate",
  NBER = TRUE,
  shading.border = FALSE,
  Enddatelab = FALSE,
  Y2lim = c(-3,33),
  Y2at = seq(0, 30, 5),
  Zeroline = TRUE
)

legend(x = "top", 
       legend = c(buscycle_legtext[4:6], "Baseline"),
       col = c(buscycle_colors[3:1], "black"),
       x.intersp = .3,
       bty = "n",
       cex = .65,
       lwd = 1.5,
       ncol = 2)
dev.off()

# Figure 17: Major asset class shares -------------------------------------

asset_classes <- c("pension", "real_estate", "corp_equ_mfs", "non_corp_eq")

for (asset_class in asset_classes) {
  
  postscript(file = file.path(figures_directory, paste0("shares_", asset_class, "_withlabs.eps")), width = 4, height = 4)
  par(mar = c(2,1,2,2))
  pinfo <- dfa_shares_stacked(plotdata = asset_shares_bynw[[asset_class]],
                              Col = buscycle_colors,
                              legend = FALSE,
                              Xlim = c(1987.5,2022.75),
                              Slab_ends = TRUE,
                              Slab_skip = c(1),
                              Y2lim = c(-10,107)
  )
  dev.off()
  
  
}

# Figure 24: Asset breakdown during recessions ----------------------------

covid_levels_stacked <- function(asset, ylim, yat) {
  postscript(file = file.path(figures_directory, paste0(asset, "_covid_levels.eps")), width = 4, height = 4)
  par(mar = c(2,1,2,2))
  dfa_levels_stacked(plotdata = covid_levels_bynw[[asset]], Col = dfa_colors[6:1], legend=FALSE, Y2lim = ylim, Y2at = yat,
                     eventdays = c(jul(20200201)), events = c("COVID crisis begins"), events.line.frac = .9)
  # NBER = TRUE)
  dev.off()
}

gr_levels_stacked <- function(asset, ylim, yat) {
  postscript(file = file.path(figures_directory, paste0(asset, "_gr_levels.eps")), width = 4, height = 4)
  par(mar = c(2,1,2,2))
  dfa_levels_stacked(plotdata = gr_levels_bynw[[asset]], Col = dfa_colors[6:1], legend=FALSE, Y2lim = ylim, Y2at = yat,
                     eventdays = c(jul(20071201), jul(20090630)), events = c("GR begins", "GR ends"), events.line.frac = c(.9, .9))
  dev.off()
}


covid_levels_stacked("bus_eq", c(-5,65), seq(0,65,10))
covid_levels_stacked("housing", c(-3,40), seq(0,40,5))
covid_levels_stacked("liquid_assets", ylim = c(-2,16), yat = seq(0,16,2))

gr_levels_stacked("liquid_assets", ylim = c(-1.5,10), yat = seq(0,10,2))
gr_levels_stacked("housing", c(-4,35), seq(0,25,5))
gr_levels_stacked("bus_eq",  c(-2.5,25), seq(0,25,5))


# Figure 25: Debt breakdown during recessions -----------------------------

covid_levels_stacked("morts", c(-1.5,12), seq(0,12,2))
covid_levels_stacked("other_liab", c(-.5,5.5), seq(0,5.5,.5))

gr_levels_stacked("morts", c(-1.5,12), seq(0,12,2))
gr_levels_stacked("other_liab", c(-.4,3.5), seq(0,3.5,.5))

# Figure 26: balance sheets during great recession - levels ---------------------------------------

for (gr in names(recession_bs)) {
  
  postscript(file = file.path(figures_directory, paste0("balance_sheet_wealth_", tolower(gr), ".eps")), width = 4, height = 4)
  # bottom, left, top, right
  par(mar = c(3, 1, 2, 2))
  pinfo <- policyPlot::rplot.bar(recession_bs[[gr]][1:5],
                                 Col = c("#4D4D4D", "#7D7D7D", "#9D9D9D", "#BDBDBD", "#CDCDCD"),#c("steelblue4", "lightskyblue2", "palegreen3", "red3" ,"indianred1"),
                                 Border = "black",
                                 Zeroline = TRUE,
                                 Y2lab = "Trillions of dollars" )
  
  # add net worth points in the center of each bar
  points(x = pinfo$x.mid, y = recession_bs[[gr]][["networth"]], pch = 20)
  
  
  dev.off()
}


# Figure 26: Balance sheets during the COVID recession - levels -----------

for (gr in names(covid_bs)) {
  
  postscript(file = file.path(figures_directory, paste0("covid_balance_sheet_wealth_", tolower(gr), ".eps")), width = 4, height = 4)
  # bottom, left, top, right
  par(mar = c(3, 1, 2, 2))
  pinfo <- policyPlot::rplot.bar(covid_bs[[gr]][1:5],
                                 Col = c("#4D4D4D", "#7D7D7D", "#9D9D9D", "#BDBDBD", "#CDCDCD"), #c("steelblue4", "lightskyblue2", "palegreen3", "red3" ,"indianred1"),
                                 Border = "black",
                                 Zeroline = TRUE,
                                 Y2lab = "Trillions of dollars" )
  
  # add net worth points in the center of each bar
  points(x = pinfo$x.mid, y = covid_bs[[gr]][["networth"]], pch = 20)
  
  dev.off()
}


# Figure 27: Leverage during recessions -----------------------------------

setcmpar()
style$legend.cex <- 1
style$freq.label.cex <- 1


dfa_buscycle_ratio_chart <- function(plot.data, Y2.lim = NULL, Y2.int = NULL) {
  par(mar = c(2,1,2,2))
  pinfo <- rplot.line(
    plotdata = plot.data,
    Col = rev(buscycle_colors),
    legend = FALSE,
    Enddatelab = FALSE,
    NBER = TRUE,
    shading.border = FALSE,
    Zeroline = TRUE,
    Y2lim = Y2.lim,
    Y2int = Y2.int
  )
}


###  Leverage


# Great recession

postscript(file = file.path(figures_directory, "leverage_gr.eps"), width = 4, height = 4)
dfa_buscycle_ratio_chart(plot.data = dfa_ratios_gr_list[["leverage"]])
dev.off()


# COVID

postscript(file = file.path(figures_directory, "leverage_covid.eps"), width = 4, height = 4)
dfa_buscycle_ratio_chart(plot.data = dfa_ratios_covid_list[["leverage"]])
dev.off()


###  Liquidity

# Great recession 

postscript(file = file.path(figures_directory, "liquidity_gr.eps"), width = 4, height = 4)
dfa_buscycle_ratio_chart(plot.data = dfa_ratios_gr_list[["liquidity"]] %>% select(-Bottom50),
                         Y2.lim = c(0,.9),
                         Y2.int = .1)
dev.off()


# Covid 

postscript(file = file.path(figures_directory, "liquidity_covid.eps"), width = 4, height = 4)
dfa_buscycle_ratio_chart(plot.data = dfa_ratios_covid_list[["liquidity"]] %>% select(-Bottom50),
                         Y2.lim = c(0,.9),
                         Y2.int = .1)
dev.off()


